/*! \file */
/*
* Copyright (c) 2008-2011, Thomas Handorf
* 
* This file is part of libScopes.
* libScopes is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
* 
* libScopes is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
* 
* You should have received a copy of the GNU Lesser General Public License
* along with libScopes.  If not, see <http://www.gnu.org/licenses/>.
*/
#include "findseeds.h"
static void aidxset(VC* a,VS* IDX, size_type i, size_type num, size_type nrc, char val);

/*!
	\brief Find the minimal seed from which the network can produce the target.
	\param net the network
	\param idx the order of the compounds, the values in this array are the indices of the compounds in the network
	\param tgt the target set of compounds as mask (original order of the network)
	\param expct number of seed compounds expected, default=0
	\param brk the last brk compounds in the idx array will not be removed, hence will be always in the seed; default=0
	\param info info struct which indicates number of scopes calculated and number of seeds found
	
	The result of this function is unique. For obtaining different seeds for a target provide different orderings in the idx array.
	Choose expct carefully. In particular underestimating expct will result in a significant performance loss.
*/
VC* findseed(Net* net, VS* idx, VC* tgt, size_type expct, size_type brk, findseed_info* info){
	size_type nsc,nsd;
	size_type i,num,nnum;
	size_type nrc=net->numC(); // number of compounds in network
	VC* sc; // takes the scopes
	VC* seed=new VC(nrc,1); // will hold the seed
	
	i=0;
   num=(expct ? (nrc-brk-i)/expct : (nrc-brk-i)/10); //set search block size according to number of expected seeds
	if (num<1) {num=1;}   //search block size should be larger than zero
	nsd=0; //number of seeds
	nsc=0; //number of scopes calculated
	
	aidxset(seed,idx,i,num,nrc-brk,0); //deleting seed compounds in current search block
	
	do {
		sc=scope(net,seed);   // get scope
		nsc++;
		if (tgt->isSubsetOf(*sc)){ //target still in scope - leave all removed seed compounds outside
			i+=num;
			// num keeps the value it has as the essential compound is in the other half search block
         num=(expct ? (nrc-brk-i)/(std::max(1,(int) (expct-nsd))) : (nrc-brk-i)/10);
         if (num<1) num=1;   //search block size should be larger than zero
         aidxset(seed,idx,i,num,nrc-brk,0); //deleting seed compounds in next search block
      } else {  //target not in seed
			if (num==1){   //single compound
				aidxset(seed,idx,i,1,nrc-brk,1);
				i++;
				nsd++;
//				printf("(%4d,%4d) ",i,nsd);
            num=(expct ? (nrc-brk-i)/(std::max(1,(int) (expct-nsd))) : (nrc-brk-i)/10);
/*				miss=expct-nsd;     // how many seed compounds do we expect in the remaining list
				if (miss<1) {miss=1;}  //should always assume that there is still one seed compound missing
				num=(nrc-brk-i)/(miss); //continue with full search block size*/
				if (num<1) {num=1;}   //search block size should be larger than zero
				aidxset(seed,idx,i,num,nrc-brk,0); //deleting seed compounds in next search block
			} else {     // still several possibilities
				nnum=num/2;   // divide search block in two - and continue checking left block
				aidxset(seed,idx,i+nnum,num-nnum,nrc-brk,1);  // return seeds of the second half of the search block
				num=nnum;
			}
		}
		delete sc;
	} while (i<nrc-brk);  //do not check the last #brk compounds - let them always in the seed
//	} while ((i<nrc-brk) && (_break==0));	
	
//	printf("scopes: %ld, seeds: %ld\n",nsc,nsd);
   if (info){
      info->numseeds=nsd;
      info->numscopes=nsc;
   }
	return seed;
}

   

/*! \brief Weighted randomization of an array
	\param IDX the index array to be randomized
	\param WGT weigths indicating the preferred sorting
	\param beta parameter of deviation from the preferred sorting
	\param num number of permutations
	\param nrc only the first nrc entries in IDX are affected
	
	Randomizes a index array while trying to keep it close to a preferred ordering. For beta=0 the index array is sorted such that weights in WGT at the positions indicated by IDX are sequencially sorted (heigher weights in the beginning of the list). For larger beta values this ordering is distorted (beta->inf results in a complete randomization). Note, the num value has to be chosen suficiently high in order to allow the randomization process to converge. This can be tested be looking at the distributions for various values of num.
*/	
void weightRand(VS* IDX, VF* WGT, double beta, size_type num, size_type nrc){
	size_type i,p1,p2,h,swap;
//	_aprint(IDX,nrc);
//	size_type nrc=IDX.size(); 
   if ((long) nrc<0) nrc=IDX->size();
	for (i=0;i<num;i++){
		p1=rand() % nrc;  // assuming RAND_MAX >> nrc !!!
		p2=rand() % nrc;
		if (p1>p2) {
			h=p1;
			p1=p2;
			p2=h;
		}  // now p1 is the lower index
		swap=0;
		if ((*WGT)[(*IDX)[p1]]>(*WGT)[(*IDX)[p2]]){
			// compounds in right order (higher weight in lower list position)
//			printf("%d %d %f\n",p1,p2,exp(((float) (*WGT)[(*IDX)[p2]]-(*WGT)[(*IDX)[p1]])/beta));
			if ((beta>0) && (rand()<RAND_MAX*exp(((float) (*WGT)[(*IDX)[p2]]-(*WGT)[(*IDX)[p1]])/beta))){
				swap=1;
			} 
		} else {
			// always swap if lower weight is on lower position -> this part always tries to keep the list sorted
			swap=1;

		}
		if (swap){
			h=(*IDX)[p1];
			(*IDX)[p1]=(*IDX)[p2];
			(*IDX)[p2]=h;  // swap compounds at corresponding list positions
		}
	}
}

/*! \brief Full randomization of an array.
	\param IDX the array to be randomized
	\param nrc only the first nrc entries in IDX are affected
*/
void fullRand(VS* IDX, size_type nrc){ //full randomization of array
	size_type i,i2,h;
   if ((long) nrc<0) nrc=IDX->size();
	for (i=0;i<nrc-1;i++){
		i2=rand() % (nrc-i);
		h=(*IDX)[i];
		(*IDX)[i]=(*IDX)[i+i2];
		(*IDX)[i+i2]=h;
	}
}

/*! \brief sorts index array according to scores
	\param idx the index array
	\param score the scores
	\param nrc only the first nrc entries in IDX are affected
	
	sorts the index array such that the indices at the begining of the index array point to large values in score array and vice versa.
*/
template <class T> void idxSort(VS* idx, array<T>* score, size_type nrc){
   if ((long) nrc<0) nrc=idx->size();
	_idxSort(idx,score, 0, nrc-1);
}
template <class T> void _idxSort(VS* idx, array<T>* score, size_type left, size_type right){
	//quicksort
	if (left>=right) return;
	size_type pvt=left+((right+1-left) / 2);
	T pivot=(*score)[(*idx)[pvt]];
	std::swap((*idx)[pvt],(*idx)[right]);
	size_type i,store=left;
	for (i=left;i<right;i++){
		if ((*score)[(*idx)[i]]>pivot){   // sorts larger elements to the beginning of the list
			if (i!=store){
				std::swap((*idx)[store],(*idx)[i]);
			}
			store++;
		}
	}
	std::swap((*idx)[store],(*idx)[right]);
	_idxSort(idx,score,left,store-1);
	_idxSort(idx,score,store+1,right);
}

/*!
	\brief array fill method for permutated array.
	\param a array to put values in
	\param IDX array containing the order of indices
	\param i start position in IDX
	\param num number of elements to use from IDX
	\param nrc maximum length of array IDX, so even if i+num is larger than nrc, only elements up to nrc-1 are used
	\param val value to be put in a
	
	puts value val in array a at indices found in array IDX between position i and i+num
*/
void aidxset(VC* a,VS* IDX, size_type i, size_type num, size_type nrc, char val){
	size_type j,max;
	
	max=i+num;
	if (max>nrc){ max=nrc; }
	
	for (j=i;j<max;j++){
		(*a)[(*IDX)[j]]=val;
	}
}
